Basics of BigDFT: first runs and managing different calculations, N2 molecule as example

This is a simple notebook that shows how to execute a simple calculation with BigDFT. You will learn how to manipulate basic DFT objects from a python script.

Atomic positions, examples of N2 Molecule

For this tutorial, we will run just a simple N2 molecule, taken from the BIgDFT database.

[1]:
from BigDFT.Systems import System
from BigDFT.Fragments import Fragment
from BigDFT.IO import XYZReader

N2 = System()
with XYZReader("N2") as ifile:
    N2["N:0"] = Fragment([next(ifile)])
    N2["N:1"] = Fragment([next(ifile)])

Every system object has a posinp representation, which is a dictionary representation of the geometry. We will use this representation later.

[2]:
print(N2.get_posinp())
{'units': 'angstroem', 'positions': [{'N': [0.0, 0.0, 0.5488], 'frag': ['N', '0']}, {'N': [0.0, 0.0, -0.5488], 'frag': ['N', '1']}]}

We now run the DFT calculation (with the default input parameters) of this molecule. This can be done by instantiating a class from the Calculators module. Then we call the run method of that calculator with the posinp of N2 as a parameter.

[3]:
from BigDFT import Calculators as C
study = C.SystemCalculator(verbose=False) #Create a calculator
log = study.run(posinp=N2.get_posinp()) #run the code

The run() method of the BigDFT.Calculators.SystemCalculator class shows in the standard output which is the command that is invoked. Then, an instance of the Logfile class is returned. Such class can be used (we will see after) to extract the information about the electronic structure of the system.

Now we can retrieve important informations on the run. See the examples below:

[4]:
log.energy #this value is in Ha
[4]:
-19.884615273244304
[5]:
log.evals[0].tolist() # the eigenvalues in Ha ([0] stands for the first K-point, here meaningless)
[5]:
[[-1.041353673208,
  -0.4926440819324,
  -0.4357815638016,
  -0.4357812404171,
  -0.3818323179544]]
[6]:
log.get_dos().plot() #the density of states
../_images/tutorials_N2_11_0.png

How to modify the Input parameters from the default one

To specify non-default input parameters to the code, we should employ a Inputfile object. One can choose a naming prefix for a run, which enables to classify the runs which are performed, and ease the subsequent postprocessing.

Methods of the class can be employed for modifying the input parameters. For instance, the XC functional can be chosen via the set_xc method. Imagine for example that you are interested in visualizing the wavefunctions output of the calculation. To do that, the suitable method of the class instance should be called. Create a new calculation set by using, for instance, the LDA prefix.

[7]:
from BigDFT import Inputfiles as I
inp = I.Inputfile()
inp.set_xc('LDA')
inp.write_orbitals_on_disk()

The wavefunctions will be present at the end of calculation, by indicating the value of the key orbitals in the output dictionary.

[8]:
#Run the code with the LDA prefix
lda = study.run(name="LDA",input=inp,posinp=N2.get_posinp())

When using a naming scheme, the output files are placed in a directory called data-LDA. In our LDA example, the wavefunctions of the system can thus be found in the data-LDA directory:

[9]:
!ls data-LDA
time-LDA.yaml                    wavefunction-k001-NR.bin.b000004
wavefunction-k001-NR.bin.b000001 wavefunction-k001-NR.bin.b000005
wavefunction-k001-NR.bin.b000002 wavefunction-rhoij.bin
wavefunction-k001-NR.bin.b000003

Here 001 means the first K-point (meaningless in this case), N stands for non spin-polarized, R for real part and the remaining number is the orbital ID. Post-processing of these files is done in the another tutorial.

In the same spirit, another calculation can be done with different parameters. Imagine we want to perform a Hartree-Fock calculation. We should simply change the XC functional used. However, we also have to specify the pseudopotential.

[10]:
inp.set_xc('HF')
inp['psppar.N']={'Pseudopotential XC': 1} #this set the default PSP for LDA
hf = study.run(name="HF",input=inp) #Run the code with the name scheme HF

Without specifying the PSP, an error will occur. The same error is specified at the end of the log file, and also in debug/bigdft-err-0.yaml:

 Additional Info:
   The pseudopotential parameter file "psppar.N" is lacking, and no registered pseudo found
   for "N"

This is because the pseudopotential is assigned by default in the code only for LDA and PBE XC approximations. Another alternative is to specify the internal PSP that might be used, taking from the default database of BigDFT. This can be done as follows:

[11]:
inp['psppar.N']={'Pseudopotential XC': 1} #here 1 stands for LDA as per the XC codes

When possible, care should be taken in choosing a pseudopotential which has been generated with the same XC approximation used. Unfortunately, at present HGH data are only available for semilocal functionals. For example, the same exercise as follows could have been done with Hybrid XC functionals, like for example PBE0 (ixc=-406). XC functionals are specified using the built in named functionals, or using the LibXC codes. Data of the calculations can be analysed.

[12]:
inp.set_xc('PBE0')
pbe0 = study.run(name="PBE0",input=inp) #Run the code with the name scheme PBE0

The variables lda, hf, and pbe0 contains all information about the calculation. This is a class Logfile which simplify considerably the extraction of parameters of the associated output file log-LDA.yaml. If we simply type:

[13]:
print(lda)
- Atom types:
  - N
- cell: Free BC
- number_of_orbitals: 5
- posinp_file: LDA.yaml
- XC_parameter: -20
- grid_spacing:
  - 0.45
  - 0.45
  - 0.45
- spin_polarization: 1
- total_magn_moment: 0
- system_charge: 0
- rmult:
  - 5.0
  - 8.0
- dipole:
  - -0.0005134858
  - -0.0005134858
  - -0.000559708
- energy: -19.884615273244282
- fermi_level: -0.3818323179544
- forcemax: 0.01156401391398
- forcemax_cv: 0.0
- gnrm_cv: 0.0001
- nat: 2
- symmetry: disabled
- No. of KS orbitals:
  - 5

we display some information about the LDA calculation. For instance, we can extract the eigenvalues of the Hamiltonian i.e. the DOS (Density of States):

[14]:
lda.evals[0].tolist()
[14]:
[[-1.041353673208,
  -0.4926440819324,
  -0.4357815638016,
  -0.4357812404171,
  -0.3818323179544]]

lda.log is the python dictionary associated to the full output file. The tutorial explains the basic methods. The following exercise (and its solution) gives some clues about it.

[15]:
print(pbe0)
- Atom types:
  - N
- cell: Free BC
- number_of_orbitals: 5
- posinp_file: PBE0.yaml
- XC_parameter: -406
- grid_spacing:
  - 0.45
  - 0.45
  - 0.45
- spin_polarization: 1
- total_magn_moment: 0
- system_charge: 0
- rmult:
  - 5.0
  - 8.0
- dipole:
  - -0.0004219248
  - -0.0004219248
  - -0.0004512606
- energy: -19.971244525934353
- fermi_level: -0.4512290316272
- forcemax: 0.05514199331393
- forcemax_cv: 0.0
- gnrm_cv: 0.0001
- nat: 2
- symmetry: disabled
- No. of KS orbitals:
  - 5

The overall structure of files in the disk

In the disk, after that the run is performed, different files are created: * input_minimal.yaml which contains the minimal set of input variables to perform the run; * log.yaml which contains the log of the run with all calculated quantities; * time.yaml and forces_posinp.yaml which we will see later.

For its I/O, bigdft uses the yaml format. If you look at the input_minimal.yaml file, you can see:

  #---------------------------------------------------------------------- Minimal input file
  #This file indicates the minimal set of input variables which has to be given to perform
  #the run. The code would produce the same output if this file is used as input.
 posinp:
   units: angstroem
   positions:
   - N: [0.0, 0.0, 0.0]
   - N: [0.0, 0.0, 1.114989995956421]
   properties:
     format: xyz
     source: posinp.xyz

In this case, only the atomic positions are indicated using a yaml format.

However, there is a one-to-one correspondence between a yaml file and a python dictionary. For this reason we will create the input parameters of this runs from dictionaries. See the useful yaml online parser page to understand the correspondence.

In the log file log.yaml, BigDFT automatically displays all the input parameters used for the calculations.

The parameters which were not explicitly given are set to a default value, except the atomic positions of course, which have to be given. As we did not specified input files here, our run is a single-point LDA calculation, without spin-polarisation.

Exercise

Compare the values of the HOMO and HOMO-1 eigenvalues for the LDA and the HF run. Change the values of the hgrid and crmult to find the converged values.

Note that, both in the LDA and in the HF calculation, a norm-conserving PSP is used.

The results can be compared to all-electron calculations, done with different basis sets, from references (units are eV) (1) S. Hamel et al. J. Electron Spectrospcopy and Related Phenomena 123 (2002) 345-363 and (2) P. Politzer, F. Abu-Awwad, Theor. Chem. Acc. (1998), 99, 83-87:

eigenvalues

LDA(1)

HF(1)

HF(2)

(Exp.)

3σg

10.36

17.25

17.31

(15.60)

1πu

11.84

16.71

17.02

(16.98)

2σu

13.41

21.25

21.08

(18.78)

The results depends, of course, on the precision chosen for the calculation, and of the presence of the pseudopotential. As it is well-known, the pseudopotential appoximation is however much less severe than the approximation induced by typical XC functionals. We might see that, even in the HF case, the presence of a LDA-based pseudopotential (of rather good quality) does not alter so much the results.

Here you can find the values from BigDFT calculation using a very good precision (hgrid=0.3, crmult=7.0). Note that 1 ha=27.21138386 eV.

eigenvalues

LDA

HF

3σg

10.40

17.32

1πu

11.75

16.62

2σu

13.52

21.30

How much do these values differ from the calculation with default parameters? Do they converge to a given value? What is the correlation for the N2 molecule in (PSP) LDA?

See here the solution.